################################################################
# Monte carlo for grouped binary data with t, t^2, t^3 terms
# Function
#
#
# Created by dave carter & curt signorino, 4/27/2007, modified by LKW on 12-2-14; modified for PA R&R on 10-6-15
################################################################

rm(list = ls())

ptim=proc.time()[3]
library(mgcv)
#library(rpvm)
library(snow)
library(foreign)
library(splines)
library(ggplot2)
library(svMisc)
library(reshape)
library(tcltk)
library(moments)


###########################################################################
################## Simulate temporal dependence ###########################
###########################################################################

gentd <- function(name, Nobs, h, b0, b1, xrange, s, sc) {
   xvar=-999
   yvar=-999
   tim=-999
   group=-999
   gg=1
   tt=1

   for (i in 1:Nobs){
      xx = runif(1, xrange[1], xrange[2])  	# uniform dist covariate
	p=1/(1+exp(-(b0+(b1*xx)+h[tt])))      	# logit link    
	yy=(runif(1)<p)                  		# y is binary
      xvar=rbind(xvar,xx)                  	# append x
      yvar=rbind(yvar,yy)                  	# append y
      tim=rbind(tim,tt)                		# append time
      group=rbind(group,gg)            		# append group
      tt=(yy==0)*(tt+1)+(yy==1)*(1)    		# update time counter
      gg=(yy==0)*(gg)+(yy==1)*(gg+1)   		# update group counter
   }

   x=matrix(xvar[2:(Nobs+1)],Nobs,1)
   y=matrix(yvar[2:(Nobs+1)],Nobs,1)
   tim=matrix(tim[2:(Nobs+1)],Nobs,1)
   group=matrix(group[2:(Nobs+1)],Nobs,1)

   x<-as.numeric(x) 
   df = data.frame(x, y, tim, tim^2, tim^3, group)     
   write.dta(df, file = paste0(name, "_", s, "_", sc, ".dta"))	
   return(df)
}

simtd <- function(name, sims, no.models, maxt, kn, sc = 1) {

   Beta.logit <- matrix(NA, sims, 2)
   Beta.td <- matrix(NA, sims, 150)
   Beta.ass <- matrix(NA, sims, 11)
   Beta.s <- matrix(NA, sims, 10)
   Beta.cp <- matrix(NA, sims, 10)

   df <- matrix(NA, sims, no.models-1)

   SE.logit <- matrix(NA, sims, 2)
   SE.td <- matrix(NA, sims, 150)
   SE.ass <- matrix(NA, sims, 11)
   SE.s <- matrix(NA, sims, 10)
   SE.cp <- matrix(NA, sims, 10)
   Hazard <- matrix(NA, 0, no.models*3)
   Sample <- matrix(NA, 0, 6)

   ptm <- proc.time()
   pb <- txtProgressBar(min = 1, max = sims, style = 3)

   for (s in 1:sims) {
      setTxtProgressBar(pb, s)	
      hazard <- matrix(NA, maxt, no.models*3)

	tddata <- gentd(paste0(name), Nobs, h, b0, b1, xrange, s, sc)
	attach(tddata)

	sample <- matrix(NA, 1, 6)
	sample[1,1:6] <- c(s, length(y), mean(y), skewness(tim), median(tim), mean(tim))

	### Hypothetical scenarios for predicted probabilities
	newd <- data.frame(cbind(tim=seq(1,maxt), tim.2=seq(1,maxt)^2, tim.3=seq(1,maxt)^3, x = round(mean(tddata$x))))
	fulld.0 <- cbind(tddata[-1], x = 0)
	fulld.1 <- cbind(tddata[-1], x = 1)

	### Flat (exponential) functional form
	logit.fit <- glm(y ~ x, family="binomial", data = tddata)
	logit.prob <- predict(logit.fit, newdata = newd, type = "response", se.fit = TRUE)
	logit.prob.0 <- predict(logit.fit, newdata = fulld.0, type = "response")
	logit.prob.1 <- predict(logit.fit, newdata = fulld.1, type = "response")
	vcv <- vcov(logit.fit)

	Beta.logit[s,1:length(logit.fit$coefficients)] <- logit.fit$coefficients
	df[s,1] <- nrow(logit.fit$model) - logit.fit$rank
	SE.logit[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:2])

	delta <- logit.prob.1 - logit.prob.0
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 
	hazard[1:length(logit.prob$fit), 1] <- logit.prob$fit
	hazard[1:length(logit.prob$se.fit), 2] <- logit.prob$se.fit
	hazard[1:length(logit.prob$fit), 3] <- apd[1:maxt,2]

	### Temporal dummies
	td.fit <- gam(y ~ s(tim, k=max(tim)) + x, data = tddata, family="binomial")
	td.xb <- predict.gam(td.fit, newd, se.fit = TRUE)
	td.prob <- 1/(1+exp(-(td.xb$fit)))
	td.prob.se <- td.prob * (1 - td.prob) * td.xb$se.fit
	vcv <- vcov(td.fit)

	td.xb.0 <- predict.gam(td.fit, fulld.0)
	td.xb.1 <- predict.gam(td.fit, fulld.1)
	delta <- 1/(1+exp(-(td.xb.1))) - 1/(1+exp(-(td.xb.0)))
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 

	Beta.td[s,1:length(td.fit$coefficients)] <- td.fit$coefficients
	df[s,2] <- nrow(td.fit$model) - length(td.fit$coefficients)
	SE.td[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:length(diag(vcv))])

	hazard[1:length(td.prob), 4] <- td.prob
	hazard[1:length(td.prob), 5] <- td.prob.se
	hazard[1:length(td.prob), 6] <- apd[1:maxt,2]

	### Cubic polynomials
	cp.fit <- glm(y ~ tim + tim.2 + tim.3 + x, family="binomial", data = tddata)
	cp.prob <- predict(cp.fit, newdata = newd, type = "response", se.fit = TRUE)
	cp.prob.0 <- predict(cp.fit, newdata = fulld.0, type = "response")
	cp.prob.1 <- predict(cp.fit, newdata = fulld.1, type = "response")
	vcv <- vcov(cp.fit)

	Beta.cp[s,1:length(cp.fit$coefficients)] <- cp.fit$coefficients
	df[s,3] <- nrow(cp.fit$model) - cp.fit$rank
	SE.cp[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:length(diag(vcv))])
	Beta.cp <- as.matrix(Beta.cp[,1:length(cp.fit$coefficients)])

	delta <- cp.prob.1 - cp.prob.0
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 
	hazard[1:length(cp.prob$fit), 7] <- cp.prob$fit
	hazard[1:length(cp.prob$se.fit), 8] <- cp.prob$se.fit
	hazard[1:length(cp.prob$fit), 9] <- apd[1:maxt,2]

	### Cubic B-splines
	s.fit <- glm(y ~ ns(tim, df = 4, knots = kn) + x - 1, family=binomial(link="logit"), data = tddata)
	newd2 <- data.frame(cbind(tim, ns(tim, df = 4, knots = kn), x = round(mean(x))))
	vcv <- vcov(s.fit)

	s.prob <- predict(s.fit, newdata = newd2, type = "response", se.fit = TRUE)
	s.prob.0 <- predict(s.fit, newdata = fulld.0, type = "response")
	s.prob.1 <- predict(s.fit, newdata = fulld.1, type = "response")

	delta <- s.prob.1 - s.prob.0
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 

	fulld <- cbind(newd2, s.prob$fit, s.prob$se.fit)
	fulld <- unique(fulld[duplicated(fulld),])
	fulld <- subset(fulld, tim <= maxt) 
	fulld <- rename(fulld, c("s.prob$fit" = "fit"))
	fulld <- rename(fulld, c("s.prob$se.fit" = "se.fit"))

	Beta.s[s,1:length(s.fit$coefficients)] <- s.fit$coefficients
	df[s,4] <- nrow(s.fit$model) - s.fit$rank
	SE.s[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:length(diag(vcv))])
	Beta.s <- as.matrix(Beta.s[,1:length(s.fit$coefficients)])

	hazard[1:length(fulld$fit), 10] <- fulld$fit
	hazard[1:length(fulld$fit), 11] <- fulld$se.fit
	hazard[1:length(fulld$fit), 12] <- apd[1:length(fulld$fit),2]

	### Automated smoothing spline
	k.no <- ifelse(max(tddata$tim)>=10, 10, max(tddata$tim))
	ass.fit <- gam(y ~ s(tim, k = k.no) + x, family="binomial", method="REML", data = tddata)
	ass.xb <- predict.gam(ass.fit, newd, se.fit = TRUE)
	ass.prob <- 1/(1+exp(-(ass.xb$fit)))
	ass.prob.se <- ass.prob * (1 - ass.prob) * ass.xb$se.fit
	vcv <- vcov(ass.fit)	

	ass.xb.0 <- predict.gam(ass.fit, fulld.0)
	ass.xb.1 <- predict.gam(ass.fit, fulld.1)
	delta <- 1/(1+exp(-(ass.xb.1))) - 1/(1+exp(-(ass.xb.0)))
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 

	Beta.ass[s,1:length(ass.fit$coefficients)] <- ass.fit$coefficients
	df[s,5] <- nrow(ass.fit$model) - length(ass.fit$coefficients)
	SE.ass[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:length(diag(vcv))])
#	Beta.ass <- as.matrix(Beta.ass[,1:length(ass.fit$coefficients)])

	hazard[1:length(ass.prob), 13] <- ass.prob
	hazard[1:length(ass.prob), 14] <- ass.prob.se
	hazard[1:length(ass.prob), 15] <- apd[1:maxt,2]
	hazard[, 17] <- s

      hazard[1:length(h[1:maxt]), 16] <- 1/(1+exp(-(b0+(b1*round(mean(tddata$x)))+h[1:maxt])))
	Hazard <- rbind(Hazard, hazard)

	Sample <- rbind(Sample, sample)

	detach(tddata)
	if (s == sims) cat("Done!\n")
   }

   close(pb)
   print(proc.time() - ptm)

   ### Combine the Betas and SEs
   Beta.td <- as.matrix(Beta.td[,1:max(length(td.fit$coefficients))])
   SE.td <- as.matrix(SE.td[,1:max(length(td.fit$coefficients))])

   Beta <- data.frame(cbind(b1, Beta.logit[,2], Beta.td[,2], Beta.cp[,5], Beta.s[,5], Beta.ass[,2], SE.logit[,2], SE.td[,2], SE.cp[,5], SE.s[,5], SE.ass[,2], df))
   names(Beta) <- c("Beta", "B.Exp", "B.TD", "B.CP", "B.S", "B.ASS", "SE.Exp", "SE.TD", "SE.CP", "SE.S", "SE.ASS", "df.Exp", "df.TD", "df.CP", "df.S", "df.ASS")
   write.dta(Beta, file = paste0(name, '_Beta_', sc, '.dta'))
	
   Sample <- data.frame(Sample)
   names(Sample) <- c("No.Sim", "N", "Y.Perc.1", "t.Skew", "t.Median", "t.Mean")
   write.dta(Sample, file = paste0(name, '_Sample_', sc, '.dta'))

   hdata <- data.frame(cbind(seq(1,maxt), Hazard))
   names(hdata) <- c("t", "EXP.h", "EXP.h.se", "EXP.apd", "TD.h", "TD.h.se", "TD.apd", "CP.h", "CP.h.se", "CP.apd", "S.h", "S.h.se", "S.apd", "ASS.h", "ASS.h.se", "ASS.apd", "true.h", "sim")
   write.dta(hdata, file = paste0(name, '_Hazard_', sc, '.dta'))
}


gentd.nph <- function(name, Nobs, h, b0, b1, xrange, s, sc, nph) {
   xvar=-999
   yvar=-999
   tim=-999
   group=-999
   gg=1
   tt=1

   for (i in 1:Nobs){

      xx = runif(1, xrange[1], xrange[2]) + tt/nph  	# uniform dist covariate
	p=1/(1+exp(-(b0+(b1*xx)+h[tt])))      	# logit link    
	yy=(runif(1)<p)                  		# y is binary
      xvar=rbind(xvar,xx)                  	# append x
      yvar=rbind(yvar,yy)                  	# append y
      tim=rbind(tim,tt)                		# append time
      group=rbind(group,gg)            		# append group
      tt=(yy==0)*(tt+1)+(yy==1)*(1)    		# update time counter
      gg=(yy==0)*(gg)+(yy==1)*(gg+1)   		# update group counter
   }

   x=matrix(xvar[2:(Nobs+1)],Nobs,1)
   y=matrix(yvar[2:(Nobs+1)],Nobs,1)
   tim=matrix(tim[2:(Nobs+1)],Nobs,1)
   group=matrix(group[2:(Nobs+1)],Nobs,1)

   x<-as.numeric(x) 
   df = data.frame(x, y, tim, tim^2, tim^3, group)     
   write.dta(df, file = paste0(name, "_", s, "_", sc, ".dta"))	
   return(df)
}

simtd.nph <- function(name, sims, no.models, maxt, kn, sc = 1, nph) {

   Beta.logit <- matrix(NA, sims, 2)
   Beta.td <- matrix(NA, sims, 150)
   Beta.ass <- matrix(NA, sims, 12)
   Beta.s <- matrix(NA, sims, 10)
   Beta.cp <- matrix(NA, sims, 10)

   df <- matrix(NA, sims, no.models-1)

   SE.logit <- matrix(NA, sims, 2)
   SE.td <- matrix(NA, sims, 150)
   SE.ass <- matrix(NA, sims, 12)
   SE.s <- matrix(NA, sims, 10)
   SE.cp <- matrix(NA, sims, 10)
   Hazard <- matrix(NA, 0, no.models*3+1)
   Sample <- matrix(NA, 0, 6)

   ptm <- proc.time()
   pb <- txtProgressBar(min = 1, max = sims, style = 3)

   for (s in 1:sims) {
      setTxtProgressBar(pb, s)	
      hazard <- matrix(NA, maxt, no.models*3+1)

	tddata <- gentd.nph(paste0(name), Nobs, h, b0, b1, xrange, s, sc, nph)
	attach(tddata)
	sample <- matrix(NA, 1, 6)
	sample[1,1:6] <- c(s, length(y), mean(y), skewness(tim), median(tim), mean(tim))

	### Hypothetical scenarios for predicted probabilities
#	newd <- data.frame(cbind(tim=seq(1,maxt), tim.2=seq(1,maxt)^2, tim.3=seq(1,maxt)^3, x = round(mean(tddata$x))))
	newd <- data.frame(cbind(tim=seq(1,maxt), tim.2=seq(1,maxt)^2, tim.3=seq(1,maxt)^3, x = 0))
	fulld.0 <- cbind(tddata[-1], x = 0)
	fulld.1 <- cbind(tddata[-1], x = 1)

	### Flat (exponential) functional form
	logit.fit <- glm(y ~ x, family="binomial", data = tddata)
	logit.prob <- predict(logit.fit, newdata = newd, type = "response", se.fit = TRUE)
	logit.prob.0 <- predict(logit.fit, newdata = fulld.0, type = "response")
	logit.prob.1 <- predict(logit.fit, newdata = fulld.1, type = "response")
	vcv <- vcov(logit.fit)

	Beta.logit[s,1:length(logit.fit$coefficients)] <- logit.fit$coefficients
	df[s,1] <- nrow(logit.fit$model) - logit.fit$rank
	SE.logit[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:2])

	delta <- logit.prob.1 - logit.prob.0
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 
	hazard[1:length(logit.prob$fit), 1] <- logit.prob$fit
	hazard[1:length(logit.prob$se.fit), 2] <- logit.prob$se.fit
	hazard[1:length(logit.prob$fit), 3] <- apd[1:maxt,2]

	### Temporal dummies
	td.fit <- gam(y ~ s(tim, k=max(tim)) + x, data = tddata, family="binomial")
	td.xb <- predict.gam(td.fit, newd, se.fit = TRUE)
	td.prob <- 1/(1+exp(-(td.xb$fit)))
	td.prob.se <- td.prob * (1 - td.prob) * td.xb$se.fit
	vcv <- vcov(td.fit)

	td.xb.0 <- predict.gam(td.fit, fulld.0)
	td.xb.1 <- predict.gam(td.fit, fulld.1)
	delta <- 1/(1+exp(-(td.xb.1))) - 1/(1+exp(-(td.xb.0)))
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 

	Beta.td[s,1:length(td.fit$coefficients)] <- td.fit$coefficients
	df[s,2] <- nrow(td.fit$model) - length(td.fit$coefficients)
	SE.td[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:length(diag(vcv))])

	hazard[1:length(td.prob), 4] <- td.prob
	hazard[1:length(td.prob), 5] <- td.prob.se
	hazard[1:length(td.prob), 6] <- apd[1:maxt,2]

	### Cubic polynomials
	cp.fit <- glm(y ~ tim + tim.2 + tim.3 + x, family="binomial", data = tddata)
	cp.prob <- predict(cp.fit, newdata = newd, type = "response", se.fit = TRUE)
	cp.prob.0 <- predict(cp.fit, newdata = fulld.0, type = "response")
	cp.prob.1 <- predict(cp.fit, newdata = fulld.1, type = "response")
	vcv <- vcov(cp.fit)

	Beta.cp[s,1:length(cp.fit$coefficients)] <- cp.fit$coefficients
	df[s,3] <- nrow(cp.fit$model) - cp.fit$rank
	SE.cp[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:length(diag(vcv))])
	Beta.cp <- as.matrix(Beta.cp[,1:length(cp.fit$coefficients)])

	delta <- cp.prob.1 - cp.prob.0
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 

	cp.fit2 <- glm(y ~ x * tim + x * tim.2 + x * tim.3, family="binomial", data = tddata)
	cp.prob2 <- predict(cp.fit2, newdata = newd, type = "response")
	
	cp.prob2.0 <- predict(cp.fit2, newdata = fulld.0, type = "response")
	cp.prob2.1 <- predict(cp.fit2, newdata = fulld.1, type = "response")
	delta2 <- cp.prob2.1 - cp.prob2.0
	apd2 <- aggregate(delta2, list(tim), sum) / aggregate(delta2, list(tim), length) 

	hazard[1:length(cp.prob$fit), 7] <- cp.prob$fit
	hazard[1:length(cp.prob$se.fit), 8] <- cp.prob$se.fit
	hazard[1:length(cp.prob$fit), 9] <- apd[1:maxt,2]
	hazard[1:maxt, 10] <- apd2[1:maxt,2]

	### Cubic B-splines
	s.fit <- glm(y ~ ns(tim, df = 4, knots = kn) + x - 1, family=binomial(link="logit"), data = tddata)
	newd2 <- data.frame(cbind(tim, ns(tim, df = 4, knots = kn), x = round(mean(x))))
	vcv <- vcov(s.fit)

	s.prob <- predict(s.fit, newdata = newd2, type = "response", se.fit = TRUE)
	s.prob.0 <- predict(s.fit, newdata = fulld.0, type = "response")
	s.prob.1 <- predict(s.fit, newdata = fulld.1, type = "response")

	delta <- s.prob.1 - s.prob.0
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 

	fulld <- cbind(newd2, s.prob$fit, s.prob$se.fit)
	fulld <- unique(fulld[duplicated(fulld),])
	fulld <- subset(fulld, tim <= maxt) 
	fulld <- rename(fulld, c("s.prob$fit" = "fit"))
	fulld <- rename(fulld, c("s.prob$se.fit" = "se.fit"))

	Beta.s[s,1:length(s.fit$coefficients)] <- s.fit$coefficients
	df[s,4] <- nrow(s.fit$model) - s.fit$rank
	SE.s[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:length(diag(vcv))])
	Beta.s <- as.matrix(Beta.s[,1:length(s.fit$coefficients)])

	hazard[1:length(fulld$fit), 11] <- fulld$fit
	hazard[1:length(fulld$fit), 12] <- fulld$se.fit
	hazard[1:length(fulld$fit), 13] <- apd[1:length(fulld$fit),2]

	### Automated smoothing spline
	k.no <- ifelse(max(tddata$tim)>=10, 10, max(tddata$tim))

	ass.fit <- gam(y ~ s(tim, k = k.no) + x, family="binomial", method="REML", data = tddata)
	ass.xb <- predict.gam(ass.fit, newd, se.fit = TRUE)
	ass.prob <- 1/(1+exp(-(ass.xb$fit)))
	ass.prob.se <- ass.prob * (1 - ass.prob) * ass.xb$se.fit
	vcv <- vcov(ass.fit)	

	ass.xb.0 <- predict.gam(ass.fit, fulld.0)
	ass.xb.1 <- predict.gam(ass.fit, fulld.1)
	delta <- 1/(1+exp(-(ass.xb.1))) - 1/(1+exp(-(ass.xb.0)))
	apd <- aggregate(delta, list(tim), sum) / aggregate(delta, list(tim), length) 

	Beta.ass[s,1:length(ass.fit$coefficients)] <- ass.fit$coefficients
	df[s,5] <- nrow(ass.fit$model) - length(ass.fit$coefficients)
	SE.ass[s,1:length(diag(vcv))] <- sqrt(diag(vcv)[1:length(diag(vcv))])
#	Beta.ass <- as.matrix(Beta.ass[,1:length(ass.fit$coefficients)])

	hazard[1:length(ass.prob), 14] <- ass.prob
	hazard[1:length(ass.prob), 15] <- ass.prob.se
	hazard[1:length(ass.prob), 16] <- apd[1:maxt,2]
	hazard[, 17] <- s

      hazard[1:length(h[1:maxt]), 17] <- 1/(1+exp(-(b0+(b1*round(mean(tddata$x)))+h[1:maxt])))
	Hazard <- rbind(Hazard, hazard)

	Sample <- rbind(Sample, sample)

	detach(tddata)
	if (s == sims) cat("Done!\n")
   }

   close(pb)
   print(proc.time() - ptm)

   ### Combine the Betas and SEs
   Beta.td <- as.matrix(Beta.td[,1:max(length(td.fit$coefficients))])
   SE.td <- as.matrix(SE.td[,1:max(length(td.fit$coefficients))])

   Beta <- data.frame(cbind(b1, Beta.logit[,2], Beta.td[,2], Beta.cp[,5], Beta.s[,5], Beta.ass[,2], SE.logit[,2], SE.td[,2], SE.cp[,5], SE.s[,5], SE.ass[,2], df))
   names(Beta) <- c("Beta", "B.Exp", "B.TD", "B.CP", "B.S", "B.ASS", "SE.Exp", "SE.TD", "SE.CP", "SE.S", "SE.ASS", "df.Exp", "df.TD", "df.CP", "df.S", "df.ASS")
   write.dta(Beta, file = paste0(name, '_Beta_', sc, '.dta'))
	
   Sample <- data.frame(Sample)
   names(Sample) <- c("No.Sim", "N", "Y.Perc.1", "t.Skew", "t.Median", "t.Mean")
   write.dta(Sample, file = paste0(name, '_Sample_', sc, '.dta'))

   hdata <- data.frame(cbind(seq(1,maxt), Hazard))
#   names(hdata) <- c("t", "EXP.h", "EXP.h.se", "EXP.apd", "TD.h", "TD.h.se", "TD.apd", "CP.h", "CP.h.se", "CP.apd", "S.h", "S.h.se", "S.apd", "ASS.h", "ASS.h.se", "ASS.apd", "true.h", "sim")
   names(hdata) <- c("t", "EXP.h", "EXP.h.se", "EXP.apd", "TD.h", "TD.h.se", "TD.apd", "CP.h", "CP.h.se", "CP.apd", "CP.NPH.apd", "S.h", "S.h.se", "S.apd", "ASS.h", "ASS.h.se", "ASS.apd", "true.h", "sim")
   write.dta(hdata, file = paste0(name, '_Hazard_', sc, '.dta'))
}


b.coverage <- function(b, b.se, true.b, level = .95, df = Inf) {

   qtile <- level + (1 - level)/2
   lower.bound <- b - qt(qtile, df = df) * b.se
   upper.bound <- b + qt(qtile, df = df) * b.se

   true.in.ci <- ifelse(true.b >= lower.bound & true.b <= upper.bound, 1, 0)
   cp <- mean(true.in.ci)
   mc.lower.bound <- cp - 1.96 * sqrt((cp * (1 - cp)) / length(b))
   mc.upper.bound <- cp + 1.96 * sqrt((cp * (1 - cp)) / length(b))

   return(list(coverage.probability = cp, MC.error.bounds = c(mc.lower.bound, mc.upper.bound)))
}

h.coverage <- function(h, h.se, t, true.h, level = .95, df = Inf) {
   qtile <- level + (1 - level)/2
   lower.bound <- h - qt(qtile, df = df) * h.se
   upper.bound <- h + qt(qtile, df = df) * h.se

   true.in.ci <- ifelse(true.h >= lower.bound & true.h <= upper.bound, 1, 0)
   cp <- aggregate(true.in.ci, list(t), mean)
   mc.lower.bound <- cp["x"] - 1.96 * sqrt((cp["x"] * (1 - cp["x"])) / length(h/maxt))
   mc.upper.bound <- cp["x"] + 1.96 * sqrt((cp["x"] * (1 - cp["x"])) / length(h/maxt))

   cover <- cbind(cp[2], mc.lower.bound, mc.upper.bound)
   colnames(cover) <- c("Coverage Probability", "Coverage: Lo", "Coverage: Hi")
   return(list(cover))
}

analtd <- function(name, b0, b1, sc = 1){
   B <- read.dta(paste0(name, "_Beta_", sc, ".dta"), convert.underscore = TRUE)
   
   perform_dat <- rep(0,11)
   dflist = c("Exp", "TD", "CP", "S", "ASS")
   for(i in 1:length(dflist)){
      beta_name <- paste0("B.", dflist[i])
      se_name <- paste0("SE.", dflist[i])
      df_name <- paste0("df.", dflist[i])

      temp_dat <- b.coverage(b = B[,beta_name], b.se = B[,se_name], true.b = B[,"Beta"], level = .95, df = B[,df_name])

      ab <- mean(abs(B[,beta_name] - B[,"Beta"]))		### Absolute bias
      mse <- mean((B[,beta_name] - B[,"Beta"])^2)		### Mean squared error
      mn.se <- mean(B[,se_name])					### Mean of 1000 SEs
      sd.beta <- sd(B[,beta_name])   				### SD of 1000 Betas

      perform_dat <- rbind(perform_dat, c(dflist[i], nrow(B), b0, b1, ab, mse, mn.se, sd.beta, temp_dat$coverage.probability, temp_dat$MC.error.bounds))
      colnames(perform_dat) <- c("Functional.Form", "N", "b.0", "b.1", "Absolute.Bias", "Mean.Squared.Error", "Mean.of.SEs", "SD.of.Betas", "Coverage", "Coverage.Lo", "Coverage.Hi")
   }
   perform_dat <- data.frame(perform_dat)
#   write.dta(perform_dat[2:6,], paste0(name, "_Performance_", sc, ".dta"))
   write.table(perform_dat[2:6,], paste0(name, "_Performance_", sc, ".csv"), sep = ",")
   return(perform_dat)
}

analh <- function(name, Nobs, b0, b1, maxt, sc = 1){
   H <- read.dta(paste0(name, "_Hazard_", sc, ".dta"), convert.underscore = TRUE)

   perform_dat <- data.frame(matrix(NA, maxt, 4))
   perform_dat[,1] <- seq(1: maxt)
   perform_dat[,2] <- nrow(H)/maxt
   perform_dat[,3] <- b0
   perform_dat[,4] <- b1

   dflist = c("EXP", "TD", "CP", "S", "ASS")
   for(i in 1:length(dflist)){
      h_name <- paste0(dflist[i], ".h")
      se_name <- paste0(dflist[i], ".h.se")

      temp_dat <- h.coverage(h = H[,h_name], h.se = H[,se_name], t = H[,"t"], true.h = H[,"true.h"], level = .95, df = Inf)

      perform_dat <- cbind(perform_dat, c(temp_dat))
   }
   perform_data <- data.frame(perform_dat)
   colnames(perform_dat) <- c("t", "No.Sims", "b.0", "b.1", "EXP.COVER", "EXP.COVER.LO", "EXP.COVER.HI", "TD.COVER", "TD.COVER.LO", "TD.COVER.HI", "CP.COVER", "CP.COVER.LO", "CP.COVER.HI", "S.COVER", "S.COVER.LO", "S.COVER.HI", "ASS.COVER", "ASS.COVER.LO", "ASS.COVER.HI")
   write.dta(perform_dat, paste0(name, "_Hazard_Coverage_", sc, ".dta"))
   return(perform_dat)
}



analapd <- function(name, Nobs, b0, b1, maxt, sc = 1){
   A <- read.dta(paste0(name, "_Hazard_", sc, ".dta"), convert.underscore = TRUE)

   perform_dat <- data.frame(matrix(NA, maxt, 4))
   perform_dat[,1] <- seq(1: maxt)
   perform_dat[,2] <- nrow(A)/maxt
   perform_dat[,3] <- b0
   perform_dat[,4] <- b1

   dflist = c("EXP", "TD", "CP", "S", "ASS")
   for(i in 1:length(dflist)){
   
	Perform_dat <- data.frame(matrix(NA, maxt, 3))
	apd_name <- paste0(dflist[i], ".apd")

	apd <- aggregate(A[,apd_name], list(A$t), mean)
	apd.ci <- data.frame(aggregate(A[,apd_name], list(A$t), na.rm = TRUE, FUN = quantile, probs = c(.025, .975)))

	Perform_dat[,1] <- apd[,2]
	Perform_dat[,2:3] <- apd.ci[,-1]

	perform_dat = cbind(perform_dat, Perform_dat)
   }
   perform_dat <- data.frame(perform_dat)
   names(perform_dat) <- c("t", "No.Sims", "b.0", "b.1", "EXP.APD", "EXP.APD.LO", "EXP.APD.HI", "TD.APD", "TD.APD.LO", "TD.APD.HI", "CP.APD", "CP.APD.LO", "CP.APD.HI", "S.APD", "S.APD.LO", "S.APD.HI", "ASS.APD", "ASS.APD.LO", "ASS.APD.HI")
   write.dta(perform_dat, paste0(name, "_APD_", sc, ".dta"))
   return(perform_dat)
}


analapd.nph <- function(name, Nobs, b0, b1, maxt, sc = 1){
   A <- read.dta(paste0(name, "_Hazard_", sc, ".dta"), convert.underscore = TRUE)

   perform_dat <- data.frame(matrix(NA, maxt, 4))
   perform_dat[,1] <- seq(1: maxt)
   perform_dat[,2] <- nrow(A)/maxt
   perform_dat[,3] <- b0
   perform_dat[,4] <- b1

   dflist = c("EXP", "TD", "CP", "CP.NPH", "S", "ASS")
   for(i in 1:length(dflist)){
   
	Perform_dat <- data.frame(matrix(NA, maxt, 3))
	apd_name <- paste0(dflist[i], ".apd")

	apd <- aggregate(A[,apd_name], list(A$t), mean)
	apd.ci <- data.frame(aggregate(A[,apd_name], list(A$t), na.rm = TRUE, FUN = quantile, probs = c(.025, .975)))

	Perform_dat[,1] <- apd[,2]
	Perform_dat[,2:3] <- apd.ci[,-1]

	perform_dat = cbind(perform_dat, Perform_dat)
   }
   perform_dat <- data.frame(perform_dat)
   names(perform_dat) <- c("t", "No.Sims", "b.0", "b.1", "EXP.APD", "EXP.APD.LO", "EXP.APD.HI", "TD.APD", "TD.APD.LO", "TD.APD.HI", "CP.APD", "CP.APD.LO", "CP.APD.HI", "CP.NPH.APD", "CP.NPH.APD.LO", "CP.NPH.APD.HI", "S.APD", "S.APD.LO", "S.APD.HI", "ASS.APD", "ASS.APD.LO", "ASS.APD.HI")
   write.dta(perform_dat, paste0(name, "_APD_", sc, ".dta"))
   return(perform_dat)
}










